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AN EVALUATION OF JET IMPINGEMENT HEAT TRANSFER CORRELATIONS 

FOR PICCOLO TUBE APPLICATION 

William B. Wright 
QSS Group, Inc. 

Cleveland, Ohio 44135 


ABSTRACT 

Impinging jets have been used for a wide variety of 
applications where high rates of heat transfer are 
desired. This report will present a review of heat 
transfer correlations that have been published. The 
correlations were then added to the LEWICE software 
to evaluate the applicability of these correlations to a 
piccolo tube anti-icing system. The results of this 
analysis were then compared quantitatively to test 
results on a representative piccolo tube system. 

NOMENCLATURE 

A Area ratio (dimensionless) 

c chord (m) 

c n spanwise distance between 

holes (m) 

c x chordwise distance between 

holes (m) 

C p heat capacity (kJ/kgK) 

d jet hole diameter (m) 

D diameter of semi-cylinder (m) 

G mass flux (kg/m 2 s) 

h Heat transfer coefficient (kW/m 2 K) 

k Thermal conductivity (kW/mK) 

m Mass flow rate per unit span (kg/ms) 

N r Number of jet rows 

Nu Nusselt number _ hd_ 

K 

Nu average Nusselt number 

q" Heat flux (kW/m 2 ) 

Pr Prandtl number _ CpVa 

K 

r radial distance from hole (m) 

Re Reynolds number _ Gd_ 

Ma 

s wrap distance (m) 

T Temperature (K) 

w slot width (m) 

z n distance from hole to wall (m) 

Correlation Constants 

a, B, b, C, c, c o , m, n, n x , n y , n z 


Greek Letters 

a impingement angle (degrees) 

/d viscosity (kg/m-s) 

p density (kg/m 3 ) 

Subscripts 

a air 

c crossflow 

j jet 

max maximum value 

o stagnation 

q heat flux 

INTRODUCTION 

The removal and/or prevention of ice on aircraft 
components is vital to aircraft performance and 
operation. Even small amounts of ice can have 
disastrous consequences. Because of this, several 
methods of ice prevention and removal have been 
designed. Methods of ice control can be arranged into 
two broad categories: anti-icing methods and de-icing 
methods. Anti-icing methods are concerned with the 
prevention or minimization of ice buildup on protected 
surfaces. De-icing methods are concerned with ice 
removal after and during ice build up. 

The primary means of preventing ice formation on 
wings and engine inlets for modern commercial 
transport aircraft is by extracting hot air from the 
compressor and blowing it on the inside surface of the 
leading edge through small holes drilled in a pipe. Ice 
accretion is prevented by supplying enough energy to 
evaporate the impinging water (evaporative anti-icing) 
or by maintaining a surface temperature above freezing 
(running wet anti-icing). This report will demonstrate 
that the amount of heat supplied can be calculated by 
correlations based upon hole and pipe diameter, the 
number of holes, the mass flow rate of the air and the 
supply temperature. 

Impinging jets are a widely used method for 
heating and cooling applications. Martin 1 and later 
Jambunathan et al. 2 provide a thorough review of heat 
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transfer data from impinging jets. Most of the research 
referenced in these works was performed using a single 
jet or a group of jets impinging on a flat plate. 
Correlations were presented in these reports for the 
average Nusselt number. While this information is 
useful, piccolo tubes impinge on the curved inside 
surface of a wing or inlet. Additionally, the effect of jet 
interaction, jet angle, and Nusselt number distribution 
are also of interest. 

Livingood and Gauntner 3 were the first to look at a 
row of jets impinging on a concave surface for a turbine 
cooling application. Similar studies were performed by 
Holl worth and Berry 4 and also Hrycak. 5 These reports 
included the first attempts at creating correlations for 
average Nusselt number. 

Gau and Chung 6 reported on the effects of both 
concave and convex surfaces and included results from 
slots as well as orifices. Sparrow and Lovell 7 studied 
the effects of jet impingement at oblique angles. 
Tawfek 8 expanded this analysis for a wider range of 
impingement angles and Reynolds numbers. Wu et al. 9 
performed a numerical study of different impingement 
angles using 3D Navier- Stokes. Goldstein and Sol 10 
looked at the effect of entrainment on heat transfer that 
can occur near the leading edge. Goldstein et al. 11 used a 
similar technique to examine variation of recovery 
factor. Florshuetz et al. 12 studied the span wise effects 
and the effect of staggered jet arrays versus inline 
arrays. 

Recently, experimental techniques have advanced 
that allow accurate measurement and correlation of 
local Nusselt number. Campbell et al. 13 and Crafton 
et al. 14 presented a methodology using temperature 
sensitive paint to obtain more accurate and more 
detailed measurements of heat transfer at different 
angles. Ichimiya 15 performed a similar study using a 
thermosensitive liquid crystal. Maurel and Solliec 16 
used laser Doppler velocimetry and particle image 
velocimetry to measure details of the jet flow field as 
well as heat transfer. Tawfek 17 used heat flux sensors to 
obtain local and average heat transfer measurements. 
He also obtained local pressure distributions using a 
Perspex disk. Huber and Viskanta 18 used the liquid 
crystal technique to examine the effect of jet spacing. 
Huang and El-Genk 19 used a heated plate to take 
measurements of local and average heat transfer. Brown 
et al. 20 used thermocouple measurements for a piccolo 
tube system but did not publish the coefficients of the 
correlation obtained. However, the form of their 
correlation is similar to the other reports which suggest 
that the other correlations may be useful. 

The report is divided into three sections. The first 
section will provide a description of the LEWICE 21 
model, with emphasis on the anti-icing physics. The 
second section will describe the specific correlations 
tested. The third section will compare the anti-icing 
results using these correlations with data taken in the 


Icing Research Tunnel (IRT) on a representative 
piccolo tube system. 

LEWICE 

The computer program LEWICE embodies an 
analytical ice accretion model that evaluates the 
thermodynamics of the freezing process that occurs 
when supercooled droplets impinge on a body. The 
atmospheric parameters of temperature, pressure, and 
velocity, and the meteorological parameters of liquid 
water content (LWC), droplet diameter, and relative 
humidity are specified and used to determine the shape 
of the ice accretion. The surface of the clean (un-iced) 
geometry is defined by segments joining a set of 
discrete body coordinates. The software consists of four 
major modules. They are 1) the flow field calculation, 

2) the particle trajectory and impingement calculation, 

3) the thermodynamic and ice growth calculation, and 

4) the modification of the current geometry by addition 
of the ice growth. 

LEWICE applies a time- stepping procedure to 
"grow" the ice accretion. Initially, the flow field and 
droplet impingement characteristics are determined for 
the clean geometry. The ice growth rate on each 
segment defining the surface is then determined by 
applying the thermodynamic model. When a time 
increment is specified, this growth rate can be 
interpreted as an ice thickness and the body coordinates 
are adjusted to account for the accreted ice. This 
procedure is repeated, beginning with the calculation of 
the flow field about the iced geometry, then continued 
until the desired icing time has been reached. 

The thermal module calculates the 2D transient 
(time-dependant) heat conduction in a body. It can 
handle multiple composite layers, where each layer can 
have different thermal properties including temperature 
dependence and anisotropy. Ice growth can be predicted 
with or without heat, and the model includes various 
routines for ice shedding and water runback. It can also 
function as an ice accretion simulator with heater power 
turned off. The original deicer program has been 
documented in previous reports. 21-23 

While this module was originally designed to 
simulate an electrothermal de-icer, a bleed air system 
can be analyzed through the use of appropriate 
boundary conditions. This method, based upon the 
model of Al-Khalil, 24 requires that the user supply either 
a heat flux or heat transfer coefficient as a boundary 
condition. For design purposes, users have indicated a 
preference for a system that instead inputs piccolo tube 
design parameters. These parameters include hole 
diameter, jet distance from the wall and jet spacing. The 
next section of this report describes empirical equations 
that have been developed which may be appropriate for 
anti-icing design. 
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LEWICE also contains a feature that assumes ID 
heat conduction in the normal direction. This 
assumption states that while the internal flow is 3D and 
the external flow and water collection are primarily 2D, 
heat conduction through the surface is primarily ID. 
Previously, the user could choose to assume 
evaporative heating or a constant (above freezing) 
surface temperature boundary in this mode. The 
calculation was performed separately from the icing 
analysis and did not affect the ice shape. It was created 
for performing quick design analysis of anti-icing 
systems. This feature has recently been extended for 
electrothermal and piccolo tube systems. The user can 
now specify the heat inputs in the same manner as the 
2D conduction model. LEWICE will then perform the 
same type of analysis (except assuming ID heat 
transfer), including the calculation of residual ice 
growth if insufficient heat is added. This mode will 
allow the user to perform more rapid preliminary design 
analysis of systems. 

PICCOLO TUBE CORRELATIONS 


The reports discussed earlier found that local and 
average Nusselt numbers were primarily functions of 
the jet Reynolds number, the jet spacing, and the 
distance to the wall. Effects of curvature and 
impingement angle were also studied. Some of the 
reports listed previously only provided a partial 
correlation. This section will present each of the full 
correlations found and their range of applicability. As 
with any correlation, care must be taken when 
extrapolating outside of the applicable range. Unless 
otherwise noted, the correlations were developed for 
arrays of circular holes impinging normally on a flat 
plate. Additionally, Reynolds number was calculated 
for circular cross-section holes by using the mass flow 
rate per unit span: 


Re = 


4 m c n 

7rN rJ u d 


z n/ d 

Q.6/4a 


Martin 1 

35 

V7 


1-2.2 VT 


1 + 0.2(z n / d, - 6)VT 


where A is the dimensionless area 



2,000 < Re < 100,000 
0.004 < A < 0.04 
2 < z n /d < 12 


Hrycak 4 5 

Flat Plate 

Nu„ = 0.763 Pr 0 39 Re°' 5 (U) 

Semi-Cylinder 

Nu a = 1 .85 Pr-^ Re 0 695 

14,000 < Re <67,000 
1.5 < z/d < 7 
D = 127 mm 


Gau and Chung 6 

Slot impinging on Concave Surface 
Re is based on the slot width, w. 


Nu = 0.25 lRe a68 | — | |^J 

wJ 


\0.15 


Nu 0 


= 0.729 Re 0 ' 5 


A 

. W 


v-0.14 



8 < Z n /w < 16 


Nu = 0.394 Re 1 


0.68 f A 


^-0.38/ Y 0 - 32 


6,000 < Re < 35,000 


8 < D c /w < 45.7 
2 < Z n /w < 8 


Tawfek 8 

Angled impingement 

S ma x represents the wrap distance measured from the 
impingement location where the maximum Nusselt 
number is located. 


Nu =0.16Re a71 pH — 



3,800 < Re < 40,000 
20° < a < 90° 

7 < z n /d < 30 
0.06 <d/D< 0.14 
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Table 1: Coefficients for inline holes 



C „ 


n v 


c 

1.18 

- 0.944 

- 0.642 

0.169 

m 

0.612 

0.059 

0.032 

- 0.02 

B 

0.437 

- 0.095 

- 0.219 

0.275 

n 

0.092 

- 0.005 

0.599 

1.04 


Table 2: Coefficients for staggered holes 




n v 

n v 

n 7 

C 

1.87 

-0.771 

-0.999 

-0.26 

m 

0.571 

0.028 

0.092 

0.039 

B 

1.03 

-0.243 

-0.307 

0.059 

n 

0.442 

0.098 

-0.003 

0.304 


2,500 < Re < 70,000 
5 < c x /d <15 (inline) 

5 < c x /d <10 (staggered) 

X 4 < c n /d < 8 

1 < z)d < 3 

Tawfek 17 

/ N-0.22/ \-0.38 

Nu = 0.453Pr^ Re a691 NM pd 

3,400 < Re <41,000 

2 < c n /d < 30 
6 < z"/d < 58 

Huber and Viskanta 18 

/ y 0.123/ y 0.725 

Nu = 0.285 Re 071 [yj 

3,400 < Re < 20,500 
0.25 < z n /d < 6 
4 < c n /d < 8 




(A (A 2 i 

r A 

3 


4 

a = 10 -4 

506 + 13.3 - -19.6 - +2.41 

- 

-0.0904 

- 



Id) \d) 1 

yd) 


yd) 


1 

( r \ f r N 

(A 

\ 3 

(A 

\ 4 

<3- 

II 

O 

k 

32-24.3 - +6.53 - -0.694 

- 

) +0.0257 

- 



\d) \d) 

\d; 


\d; 

1 



).0904 





BLEED AIR RESULTS 

The model used for comparison in this report is 
shown in Figure 1 . A two week entry was performed in 
the Icing Research Tunnel (IRT) in Sept. 2003 by Lee, 
Addy and Broeren 25 for the purpose of documenting 
residual ice shapes from a piccolo tube system. It is not 
the intent of the paper to report on the full findings of 
that test entry. This report will only present preliminary 
thermal data for comparison with the correlation model. 



Figure 1: Piccolo Tube Test Model 


The model had an average chord length of 1.52 m 
with a slight 8° sweep. The airfoil profile is similar to 
the NACA 23014 model used in previous IRT entries as 
shown in Figure 2. Thermocouples were placed at two 
different spanwise sections on the inside surface of the 
airfoil at the locations. A schematic of the 
thermocouple locations relative to the hole locations is 
shown in Figure 3. 



Figure 2: Profile of Airfoil Used in Comparison 

Section F TCs Section G TCs 

' TL5 * TL15 

- TL2 TUB- 


In Figure 3 above, the spanwise direction is 
horizontal while the leading edge has been unwrapped 
to show the thermocouple locations. The locations 
where the jets impinge the surface are shown as open 
circles. The holes are 1.32 mm in diameter and the 
distance between holes is 6.6 cm or 50 hole diameters. 
In the comparisons below, one complete set of jets 
spanwise was simulated. The center row of holes 
resides on the airfoil leading edge and in the 
comparisons will be referred to as Jet 1. The lower 
surface jet row was labeled Jet 2 and the upper surface 
jet row was labeled Jet 3. 

The spanwise distance between the two rows of 
thermocouples is greater than that shown in the diagram 
above, although the relative distance from a 
thermocouple to a hole is correctly displayed. The 
diagram was presented in this fashion so that both 
instrumentation rows could be shown. The 
thermocouple locations were selected such that the 
piccolo flow did not impinge directly. The 
thermocouple coordinates relative to the impingement 
location of the leading edge hole is shown in Table 3. 


Tab 

e 3: Thermocouple Locations 

TC# 

x(m) 

y(m) 

z(m) 

TL5 

0.0483 

-0.0390 

-0.0058 

TL15 

0.0467 

-0.0385 

-0.0038 

TL2 

0.0272 

-0.0316 

-0.0287 

TL18 

0.0259 

-0.0311 

0.0193 

TL4 

0.0083 

-0.0200 

-0.0135 

TL16 

0.0075 

-0.0191 

0.0041 

TL17 

0.0000 

0.0005 

0.0112 

TL3 

0.0000 

0.0008 

-0.0211 

TL22 

0.0033 

0.0169 

0.0112 

TL9 

0.0035 

0.0175 

-0.0211 

TL23 

0.0105 

0.0318 

0.0193 

TL10 

0.0112 

0.0330 

-0.0135 

TL24 

0.0204 

0.0453 

0.0041 

TL8 

0.0216 

0.0467 

-0.0287 

TL21 

0.0321 

0.0576 

0.0269 

TL11 

0.0336 

0.0588 

-0.0058 


TL4 


TL16 - 


TL3 


- TL9 
• TL10 
TL8 

■ TL11 


^ TL17 
TL22 ■ 3 *+ ■ 
TL23 • 

- TL21 
TL24 * 


Figure 3: Thermocouple Locations on 
Piccolo Tube Model 


The de-icing system also contained an inner liner 
to increase the heat transfer downstream of the tube. A 
cross-section view of the airfoil with this liner is shown 
in Figure 4. In LEWICE, the liner was modeled by 
increasing the flow velocity using the Bernoulli 
equation. The increased velocity was used in the 
correlation to get the heat transfer in this region. It 
should be noted that even if the Bernoulli assumption is 
correct, the modeling of the liner in this manner is 
wrong, since the correlation was developed for 
impinging jets and not for crossflow situations. 
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Table 5: Piccolo Inputs 
for Correlations 


Re 

11870 

z n /d (Jet 1) 

6.9 

z n /d (Jet 2) 

5.8 

z n /d (Jet 3) 

8.1 

c/d 

D 

50 



Figure 4: Cross-section of Piccolo tube with Liner 

Two cases from the IRT entry were selected to 
demonstrate the capabilities of the correlation. A warm 
temperature condition and a cold temperature condition 
with the same bleed air flow rate were selected from the 
test matrix. The conditions are shown in Table 4. The 
mass flow rate given for these conditions represented 
typical values for this test. 


Table 4: Test Conditions 



Case 1 

Case2 

T(°F) 

23.1 

-19.1 

V (kts) 

115 

110 

MVD 

29 

29 

LWC 

0.87 

0.67 

AOA 

3 

3 

time (min) 

22.5 

8.6 

m„, r (lb/ft-s) 

0.01 

0.01 

Tmi (°F) 

350 

350 


Using the mass flow rate and the piccolo tube 
geometry, values that can be input into the various 
correlations can be calculated. Table 5 summarizes 
these inputs. The values in this table were used in the 
Goldstein 11 correlation to calculate Nusselt number. The 
heat transfer coefficients from this calculation are 
shown in Figure 5 at three span wise values: over the 
centerline of Jet 1, over the centerline of Jets 2 and 3, 
and at the span location between the jets. The maximum 
heat transfer coefficient calculated was just over 
1 kW/m 2 K (176 BTU/ft 2 hr°F) and decreases rapidly 
away from the jet centerline. In the region between jets, 
the maximum heat transfer coefficient was under 
0.3kW/m 2 K which demonstrates that significant cooling 
occurs span wise. The heat transfer was increased in the 
liner region as shown by the jump in the curves. 


1.2 



-020 - 0.15 - 0.10 - 0.05 0.00 0.05 0.10 0.15 0.20 

-0.2 J s/c 


Figure 5: Heat Transfer Coefficients at Selected 
Spanwise Locations 

Calculations of surface temperature were then 
performed at each of the spanwise locations where a 
thermocouple measurement was recorded. Several other 
spanwise locations, ranging from the centerline of Jet 1 
to the centerline of Jets 2 and 3, were also calculated to 
determine the maximum surface temperature and to 
estimate the temperature gradient. A total of 100 
spanwise cases were calculated for each condition. 
Figure 6 shows the temperature distribution at the three 
spanwise locations where the heat transfer coefficients 
were plotted. 



Figure 6: Temperature Variation for Casel 


NASA/CR — 2004-2 12917 


6 


The limitations of the ID conduction 
approximation are apparent in this figure. For example, 
the surface temperature for the Jet 1 centerline location 
drops over 90°C (from over 350 K to 270 K) as soon as 
the simulation evaporates all of runback water. The 
increase in surface temperature over the liner is 
considered to be an artifact of the model. The calculated 
temperature gradients in both the chordwise and 
spanwise directions appear larger than expected. 

The spanwise variation in temperature also creates 
a difference in the calculated residual ice shape, shown 
in Figure 7. 



Figure 7: Residual Ice Shape for Casel 


The ridge ice location is predicted well by the 

model, however the mass of residual ice is 

underpredicted. There are several factors that could 
have caused this result. Primarily, the piccolo heat 
transfer coefficient could be overpredicted by the 

correlation. Additionally, it is acknowledged that the 

increased heat transfer from the liner would not be 
modeled by these correlations. 

The results from each of the 100 spanwise cuts can 
be merged into a single file for display purposes. The 
results of this process, taking advantage of the periodic 
spanwise boundary condition, is shown in Figures 8 
and 9. 


Figure 9: Lower Surface Ice for Casel 


Finally, the predicted temperatures at each 
thermocouple can be extracted from the individual runs 
and plotted against the experimental data. This result is 
shown in Figure 10. In this plot, the temperatures are 
plotted solely as a function of wrap distance even 
though all of the thermocouples are not at the same 
spanwise locations. The results show that LEWICE 
predicts higher temperatures than the experimental data. 
This result is consistent with the residual ice shape 
result. Higher surface temperatures will lead to higher 
amounts of evaporation and a smaller runback shape. It 
is estimated that the correlation would need to be 20% 
lower in order to match the experimental temperature 
profile. 



-0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.05 

s/c 


Figure 10: Temperature Comparison for Casel 


The same factors that affected the first case can also be 
seen in the results from Case 2. Figure 11 shows the 
surface temperature distribution for this condition. This 
case shows a larger variation in surface temperature due 
to the liner. This result suggests that the effects of the 
Figure 8: Upper Surface Ice for Casel liner were overpredicted by this model. 
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The differences and percent differences were 
calculated by the following equations: 

abs.diff =1 T exp — T LEWICE I 

IT -T I 

%dijf = — ^ LEWICE 

^exp,max ^oo 

For Case 1 , the average temperature difference was 
5.1 °C (9.2 °F) or 18% of the maximum temperature 
rise. 



Figure 11: Surface Temperature Distributions for 
Case 2 

This result also shows that the temperature between 
the jets and the liner is only slightly above freezing. A 
slightly lower piccolo heat transfer coefficient could 
result in ice formation at this location. It should also be 
noted that LEWICE was modified for these runs such 
that the external boundary layer was kept laminar over 
the heated region. If this assumption does not hold, a 
higher external heat transfer coefficient can result. 
These factors help to explain the residual ice shape 
prediction in Figure 12. 



Figure 12: Residual Ice Shape for Case 2 


In this case, LEWICE predicts an aerodynamic ally 
severe ice ridge similar to but larger than the 
experimental ice shape. LEWICE places the ice ridge 
behind the liner boundary rather than in front as shown 


by the experimental data. This result shows that an 
improved model for the liner can result in more 
accurate predictions. The span wise cuts were again 
merged to produce a pseudo-3D view of the LEWICE 
results. The upper and lower surface ice accretions are 
shown in Figures 13 and 14, respectively. 



Figure 13: Upper Surface Ice for Case2 



Figure 14: Lower Surface Ice for Case2 


The predicted temperatures at each thermocouple 
location were again extracted from the individual runs 
and plotted in Figure 15. In this case, the LEWICE 
predictions appear to be closer to the measured 
experimental data. As such, the amount of residual ice 
is closer to the experimental data. However, since the 
model used for the liner increased temperatures 
downstream, the ice did not form until after the runback 
water had passed the liner. 
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320 

315 



28C|'P T T T T T , 

-0.06 -0.04 -0.02 0 0.02 0.04 0.06 

s/c 

Figure 15: Temperature Comparison for Case 2 

In this case, the average temperature difference was 
4.3°C (7.8°F) while the percentage difference was 
10.5%. This percentage is much lower due to the larger 
denominator used. 

CONCLUSIONS 

A literature search was performed which yielded 
several candidate correlations that could be applicable 
to piccolo tube anti-icing. Representative values for 
heat transfer coefficient were calculated from one such 
correlation and applied to a piccolo tube system for 
which temperature data exists. Surface temperatures 
and residual ice shape were compared to the 
experimental test data for two cases. The results 
showed that the correlation overpredicted surface 
temperature, which resulted in different ice residual 
patterns from the experimental shapes. However, the 
technique shows promise for rapidly obtaining first 
order estimates of piccolo tube performance. A separate 
correlation will be needed for handling the inner liner. 
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